Stochastic effects in a seasonally forced epidemic model 
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The interplay of seasonality, the system's nonlinearities and intrinsic stochasticity is studied for 
a seasonally forced susceptible-exposed-infective-recovered stochastic model. The model is explored 
in the parameter region that corresponds to childhood infectious diseases such as measles. The 
power spectrum of the stochastic fluctuations around the attractors of the deterministic system that 
describes the model in the thermodynamic limit is computed analytically and validated by stochastic 
simulations for large system sizes. Size effects are studied through additional simulations. Other 
effects such as switching between coexisting attractors induced by stochasticity often mentioned in 
the literature as playing an important role in the dynamics of childhood infectious diseases are also 
investigated. The main conclusion is that stochastic amplification, rather than these effects, is the 
key ingredient to understand the observed incidence patterns. 
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I. INTRODUCTION 

Epidemic spread in human populations is a complex 
phenomenon whose comprehensive modeling has been a 
standing challenge for many years [lj . Several ingredients 
such as host population contact structure, host hetero- 
geneity, transmission mechanisms, interplay between im- 
mune response and pathogen evolution or demographic 
and environmental factors have been identified as play- 
ing an important role in the short and in the long term 
behavior of infection spread Scarcity of data, es- 

pecially for long term behavior, and model parameters 
that are hard to measure accurately add to the complex- 
ity of the problem, making it difficult to identify the key 
ingredients of a parsimonious model for a given disease 
or class of diseases [13, EH • 

Childhood infectious diseases have often been taken as 
a case study and model testing ground, because decades 
long of fairly well time resolved data records are avail- 
able, on one hand, and because of their different phe- 
nomenology despite the similarities in contagion mecha- 
nisms and in infectious, latent and immunity waning typ- 
ical times [16| . The common modeling approach for this 
class of diseases is a SEIR (susceptible-exposed-infective- 
recovered) compartmental model, with a periodic forcing 
that represents seasonal environmental or social factors 
that influence the transmission of the disease lil- 20j|. 
Deterministic models based on this approach, where the 
role of stochasticity is merely that of making the system 
switch between coexisting attractors, successfully repro- 
duce the main features of the observed time series for 
measles incidence, but fail conspicuously to model the 
behavior of other childhood diseases that exhibit non- 
seasonal sustained oscillations on several long term data 
records [2]], [22| ■ This failure has been addressed in the 
literature by claiming that this different dynamics may 
be either the result of more realistic latent and infectious 



period distributions or the evidence of stochastic effects 
that would show up as noisy oscillations with the same 
frequency as the damped oscillations of the determinis- 
tic system 15|, l22J ■ The idea that stochastic effects may 



play a more fundamental role has also been discussed in 
the epidemiological literature [14J, l23[ , prompting several 
recent analytical studies, all of which deal with unforced 
systems [Ij-Sj. 



Sustained oscillations typical of the incidence patterns 
of childhood infectious diseases are one of the features 
of the long term behavior of these unforced stochastic 
models. The power spectrum of the fluctuations around 
the endemic equilibrium computed analytically using van 
Kampen's system size expansion has well defined res- 
onant like peaks [24]], which means that for moderate 
system sizes demographic stochasticity will generate sus- 
tained noisy oscillations, a phenomenon dubbed stochas- 
tic amplification when it was first studied in ecological 
and epidemiological models [24, 27 1. Indeed, it has been 
shown that intrinsic noise enhanced by correlations in un- 
forced epidemic models may give rise to oscillations that 
are comparable to those due to seasonal forcing in de- 
terministic systems 28|-3(J- It has also been shown that 
the dominant frequency of these stochastic fluctuations 
may differ significantly from the characteristic frequency 
of the deterministic system, especially in the presence 
of correlations 0, 0, [3l| . Adding seasonality to these 
unforced systems may also give rise to significant non 
trivial effects in their long term behavior. Therefore, 
the analytical results for the fluctuations power spec- 
trum must be extended to the corresponding periodically 
forced models in order to assess the role of stochastic ef- 
fects in childhood diseases epidemiology through the in- 
terplay between seasonality, the system's nonlinearities 
and intrinsic stochasticity. 

The method developed in Ref. [27| was extended to 
fluctuations around cycles of forced or unforced systems 



2 



in a series of recent papers 32|, |33J. Here, we apply this 
method to a seasonally forced SEIR model in a realistic 
parameter region for childhood infectious diseases where 
different attractors exist or coexist. Since very little is 
known about the amplitude of seasonal forcing, we leave 
this as a free parameter in our study and consider sepa- 
rately the low, intermediate and high forcing regimes. In 
all cases, we find an excellent agreement between the an- 
alytical power spectra and the results of stochastic simu- 
lations. We use the latter to assess the role of competing 
attractors in explaining the observed time series of child- 
hood diseases epidemics, and we use the former to predict 
the number and position of the dominant non-seasonal 
peaks as a function of the epidemiological parameters. 



II. THE SEASONALLY FORCED 
DETERMINISTIC SEIR EPIDEMIC MODEL 

In this section, we will review a seasonally forced de- 
terministic SEIR model and its dynamics. Following 
most of the literature on childhood infectious diseases' 
modeling, we will take measles as a paradigmatic exam- 
ple throughout the paper, and the behavior of the sys- 
tem will be illustrated in the relevant parameter region 
3E3, EH- Extensions of the SIR/SEIR model have also 
been considered, especially in the mathematical litera- 
ture, in connection with acute infectious diseases' mod- 
eling. These extensions may include, for instance, higher 

immi- 
How- 



order nonlinearities in the infection term 34 36 



gration of infectives [2J, |37| or age structure 
ever, it is generally accepted that the seasonally forced 
SEIR model should capture the main features of the dy- 
namics of childhood infections. 

Consider then a homogeneously mixed population of 
constant size consisting of four classes of individuals: 
susceptibles (healthy individuals capable of contracting 
the infection), exposed (infected but not yet infectious 
individuals), infectives (infectious individuals capable of 
transmitting the infection) and recovered (permanently 
immune individuals). The dynamics of the SEIR model 
is governed by the following processes: 

1) The susceptible individuals catch the infection from 
the infective individuals at a time dependent contact (or 
transmission) rate Pit). To incorporate seasonality in 
transmission of the infection, the contact rate is assumed 
to be a periodic sinusoidal function with period 1 year 



/3(i)=A)(l + /3i cos27rf), 



(1) 



where t is measured in years. In Eq. (p}, ft <E [0,1] is 
the amplitude of the seasonal variation in transmission 
and /?o > is the time-averaged contact rate. This form 
of the seasonal contact rate captures the first mode of 
periodic change in disease transmission due to the school 
terms and holidays 39]. 



2) The exposed individuals become infectious at the 
rate \- Hence, 1/x is the average latent period of the 
disease. 

3) The infective individuals recover at the rate 7 be- 
coming permanently immune. The average infectious pe- 
riod of the disease is 1/7. 

4) All individuals are subjected to the per capita death 
rate /x which is equal to the birth rate of the population. 
The average lifespan is given by l//i. 

The seasonally forced deterministic SEIR model can 
now be written as follows 



ds 

~dt 

de 

~dt 

di 
di 



M(l - s(t)) - A)(l + ft cos2Trt)s(t)i(t), (2) 
ft(l + Pi cos 2irt)s(t)i{t) - ( X + n)e{t), (3) 
X e{t) - (7 + (4) 



where s(t), e(t) and i(t) denote the fraction of suscepti- 
ble, exposed and infective individuals, respectively. Note 
that the equation for the fraction of recovered individu- 
als, r(t), is redundant since s(t) + e(t) + i(t) + r(t) = 1, 
where the size of the total population was normalized to 
unity. 

The unforced deterministic SEIR model is recovered 
by setting ft = [JEl: 



ds 

di 
de 

~di 

di 
dt 



n(l-s(t))-p s(t)i(t), 
p s{t)i{t)-{ X + »)e(t), 
X e(t) - (7 + l*)i(t). 



(5) 
(6) 
(7) 



Both model ©-(ED) an d model (Et-Gl) have been ex- 
tensively investigated in the literature Hi 20 . Here, we 
will review well known facts about these models that are 
relevant for the present study. The dynamics of the un- 
forced deterministic SEIR model (JSJ-Q depends on the 
basic reproductive ratio of the infection [l|, llll 



i?n — 



PoX 



(x + /-i)(7 + M)' 



(8) 



defined as the average number of secondary cases caused 
by an infectious individual in a totally susceptible pop- 
ulation in one infectious period. The stability analysis 
shows that Eqs. ©-(JT]) have an asymptotically stable 
endemic equilibrium 



{s ,e ,1 ) = — , , 

\Ro PoX Po 

provided Rq > 1. The trivial steady state 



(9) 



(s*,e*,i*) = (1,0,0) (10) 
is stable for Rq < 1 and unstable for Rq > 1. 
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Fixed parameters 


Rate of disease onset 


X 


35.84 year" 1 


Average latent period 


Vx 


10.18 days 


Rate of recovery 


7 


100 year" 1 


Average infectious period 


1/7 


3.65 days 


Birth/death rate 


/" 


0.02 year" 1 


Average lifespan 


i/n 


50 years 


Average contact rate 




1575 year" 1 


Basic reproductive ratio 


Ro 


15.74 


Variable parameter 


Amplitude of seasonal forcing 


Pi 


0.02 (T = 1) 
0.05 (T = 1) 
0.10 (T = 1,3) 
0.12 (T = 2,3) 
0.2 (T = 2) 



TABLE I: Parameter values for measles that will be used in 
this paper. According to Eq. ([8} this set corresponds to 
Ro = 15.74. The amplitude of seasonal forcing will be varied 
so that the solutions of Eqs. ©-Q exhibit stable limit cycles 
of period T indicated in the parenthesis next to the /3i value. 
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If the contact rate varies seasonally according to Eq. 
(UJ) then a wide range of dynamical behavior of the sea- 
sonally forced deterministic SEIR model ©-(31) is possi- 
ble, depending on Rq and on the amplitude of seasonal 
forcing j3\ . In what follows, we will restrict our analysis to 
the case of measles and use the parameter values given in 
Refs. [13, El. The birth/death rate will be fixed at 0.02 
year -1 which corresponds to the average lifespan of 50 
years. The average latent and infectious periods will be 
equal to 10.18 days and 3.65 days, respectively, yielding 
13.83 days for the average time interval between infection 
and recovery. These values correspond to 35.84 year -1 
for the disease onset rate and 100 year -1 for the recovery 
rate. Finally, the average contact rate is adjusted to be 
1575 year -1 so that the basic reproductive ratio is 15.74. 
The remaining parameter, the forcing amplitude, will be 
given different values in the interval [0.02, 0.2] (estimates 



1J|). The summary of the 



for /3i can be found in Ref. 
parameter values is shown in Table fl] 

The bifurcation analysis of the seasonally forced deter- 
ministic SEIR model ©-© for the fixed set of parameter 
values considered above and Bi as a single free param- 
eter can be found in Refs. [It], [l9|]. In Ref. [l8| an 
analysis of the same model was performed for variable 
Pi and Rq — 18 (corresponding to f3 — 1800 year" 1 ). 
For the interested reader, we recommend Ref. [20[ where 
more general bifurcation diagrams of Eqs. ©-© were 
computed with two free parameters Rq and Pi and the 
remaining parameters held constant as in Table HI 

For the parameter values of Table Q] and variable Pi , 
a brief summary of the bifurcation diagram is as follows 
17, lj|. If (3\ is positive but small, a stable limit cycle of 



FIG. 1: The fraction of infective individuals as a function of 
time for periodic solutions of Eqs. ©-(fj)). 



period 1 bifurcates from the endemic equilibrium point 
©. As Pi is increased monotonically, it is found that at 
a value of P[ « 0.11479 a stable branch of period 2 bifur- 
cates from the period 1 branch, and the period 1 branch 
becomes unstable for Pi > P[. Additionally, in some 
range of Pi these period 1 and period 2 limit cycles coex- 
ist with a pair of limit cycles of period 3 (one stable and 
one unstable) appearing from a saddle-node bifurcation. 
Finally, in a very narrow window of Pi the period 2 and 
period 3 branches exhibit a cascade of period-doubling bi- 
furcations as pi increases, leading to chaotic epidemics. 
The full bifurcation diagram contains yet another stable 
attractor of period 4 that coexists with stable period 2 
and period 3 limit cycles, however this attractor has a 
very small basin of attraction and it is hard to spot both 
in numerical integrations of Eqs. ©-© and in stochastic 
simulations. 

Thus, in the realistic parameter region there are three 
main stable attractors, namely stable limit cycles of pe- 
riod 1, 2 and 3. We will consider typical values of 
Pi for which i) only a limit cycle of period 1 exists 
(pi = 0.02,0.05); ii) limit cycles of period 1 and 3 co- 
exist (Pi — 0.1); iii) limit cycles of period 2 and 3 coexist 
(Pi = 0.12); iv) a limit cycle of period 2 exists (Pi = 0.2). 
From now on the reader can refer to Table U where the 
periods T of the (co)existing limit cycles are indicated 
for each value of Pi. Note that the regimes ii) and iii) 
lie in the vicinity of the bifurcation point where a stable 
limit cycle of period 2 is born (P[ w 0.11479). 
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The fraction of infectives for stable periodic solutions 
of the seasonally forced SEIR model ©-Q is shown in 
Fig. [T] The solid curves were computed by numerical in- 
tegration of Eqs. ©-Q for the parameters values given 
in Table U and three values of the forcing amplitude. In 
all cases, initial conditions were chosen on a cycle so as to 
avoid transients. The annual epidemic cycle [/3\ = 0.05, 
Fig. Q] (a)] corresponds to a small amplitude stable limit 
cycle of period 1 bifurcating from the endemic equilib- 
rium ((9]) of the unforced SEIR model ©-O). This cycle 
is present in the system if < /3i < 0.11479. The equilib- 
rium value of infectives, i* , in the unforced model ©-([7]) 
is also shown (the dashed line). The biennial epidemic 
cycle representing the alternating years of high and low 
incidence can be found for 0.11479 < /3% < 0.2 (here we 
do not consider values of fi\ > 0.2 where a stable limit 
cycle of period 2 is still present because such levels of 
forcing are considered unrealistically high). A plot of 
typical biennial epidemic is shown for f3\ — 0.2 in Fig. 
[I] (b). Such a behavior of measles is in accordance with 
the data from the New York City and from England and 
Wales recorded between 1950 and 1968 (before vaccina- 
tion) [3 Il6l . |38j. As for other stable attractors coex- 
isting with stable limit cycles of period 1 and 2 these 
are, mainly, limit cycles of period 3 as mentioned in this 
section before. A typical limit cycle of period 3 is char- 
acterized by higher amplitude outbreaks and by very low 
minima of infective incidence, see Fig. [T] (c) where the 
behavior of the solution for the fraction of infectives is 
shown for /3i =0.1. 



III. THE SEASONALLY FORCED STOCHASTIC 
EPIDEMIC MODEL 

In the seasonally forced stochastic SEIR model, a dis- 
crete population of the constant size N is divided into the 
classes of susceptibles (S), exposed (E), infectives (I) and 
recovered (R). At any instant of time, the total number 
of individuals equals N thus only three of the classes are 
independent. Denoting numbers of susceptible, exposed 
and infective individuals by mi, mi and m 3 , respectively, 
and considering the processes postulated in Section II, we 
can obtain the transition rates for the stochastic model 
as follows. 

1) The infection process takes place between a suscep- 
tible and an infective individual at the time dependent 
contact rate f3(t) and results in the transition of the sus- 
ceptible to the exposed class, S + I —I E + I. The 
corresponding transition rate is 



T-Tni-l,m2+l,m 3 _ 
'7771,7772,7773 



TTl~\ 

/3 (l + /3 1 cos27rt)-^m 3 , (11) 



where the superscript and the subscript of T denote the 
final and the initial states of the system. 



2) The transition of an individual from the exposed 
class to the infective class, E — ^> /, occurs at the rate \- 
The transition rate of this process equals 

-mi ,7712 — 1,7713+1 



V 

' TTli ,771.2 ,"1-3 



Xm 2 - 



(12) 



3) The transition rate of the recovery of an infective 
individual occurring at the rate 7, / — R, is given by 



>-j-m 1 ,m2,m 3 - 
' mi ,m2 ,7773 



7TO3. 



(13) 



The recovered individuals are permanently immune. 

4) Finally, there are four transition rates associated 
with the linked birth, R — > S, and death processes, 



S 



R, E — > R and / — > R, occurring at the rate /x: 



T-m 1 +l,m ! ,ra 3 
'mi ,m2 ,7713 

T-mi—1, m 2 ,m 3 

' 7771,7772,7773 
/ 7~777i ^7772 — 1,771 3 
' 7771,7772,7773 



r 



- f m 1 ,m 2 ,m 3 - 
mi ,7Ti2 ,m3 



(imi, 
(j,m2, 
/jm 3 . 



(14) 
(15) 
(16) 
(17) 



Note that the initial and final states for the recovery of 
an infective individual and the death of an infective in- 
dividual are equal, see Eqs. (fl~3|) and (fl~7|) but the cor- 
responding transition rates are different. To distinguish 
between the latter we use a prime in Eq. (|17[) . 

The dynamics of the stochastic system determined by 
Eqs. (fl"l"1) - ([TTl) is completely described by the master 
equation [4Ct |41| . Given the initial and boundary condi- 
tions, this equation expresses the evolution of the prob- 
ability of having a system in state with mi susceptibles, 
m,2 exposed and m 3 infectives, V(rn,t), at any positive 
time. Denoting m as the shorthand for three numbers 
mi, ma, m 3 this equation is written in the following form: 



dP(m,t) 
dt 



£ \pV(rn',t)-TZ'v{m,t)\ (18) 



m'^m 



The positive term on the r.h.s. of Eq. (|T5)) is the increase 
of V(m,t) due to the transitions from states m' to state 
m, while the negative term is the decrease due to the 
transitions from state m to states m! . The term with m 
equal to m! is excluded because its contribution is zero. 

Substituting Eqs. ([IT ]) -([T7 jl into Eq. JISJ one obtains 
the master equation for the seasonally forced stochas- 
tic SEIR model (|46|) (see the Appendix). This equation 
is in fact a system of mi • ■ m 3 coupled differential- 
difference equations similar to those used in studies of 
simple birth-death stochastic processes [H, 42|. In that 
case, an exact solution can be found in terms of generat- 
ing functions or by other means but Eq. (l46l) with three 
variables mi , mi and m 3 is too complicated to be solved 
exactly. To get insight into the dynamics of the season- 
ally forced stochastic SEIR model for large system sizes, 
we will apply a method due to van Kampen (4fj| . Sim- 
ilarly to other studies 32, 33|, we will expand Eq. (|4"6")l 
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in powers of the inverse of the system size N. First, let where 



us rewrite Eq. (|46|) in terms of the step operators ej , 
where j — 1,2,3, defined by their action on a smooth 
function /(mi, m 2 , m 3 , t) [401 ] : 

ef 1 f(mi,m 2 ,m 3 ,t) = f(m 1 ±l,m 2l m 3 ,t), (19) 



+1 



/ (m 1 ,m 2 ,m 3 ,t) = /(mi,m 2 ± l,m 3 ,t), (20) 
^f 1 f(m 1 ,m 2 ,m 3 ,t) = f(m 1 ,m 2 ,m 3 ±l,t). (21) 



Using Eqs. (fT^ - (f2Tj) the master Eq. (|46|) transforms into 
Eq. (|j7"|) given in the Appendix. 

For the seasonally forced stochastic SEIR model the 
expansion is made around a deterministic periodic solu- 
tion of Eqs. d21)-(|4|) which is a stable limit cycle of period 
T (see Ref. (33|). Thus, we make a transformation from 
the discrete variables rrij(t) to the continuous variables 
Xj(t), where j = 1, 2, 3, according to the equations 



m\(t) 
m 2 (t) 
m 3 (t) 



Ns(t) 
Ne(t) 
Ni(t) ■ 



Nxi(t), 
Nx 2 (t), 
Nx 3 (t), 



(22) 
(23) 
(24) 



where s(t),e(t) and i(t) denote the deterministic trajec- 
tory of Eqs. ©-© and xi(t),x 2 (t) and x 3 (t) are the 
corresponding stochastic fluctuations about it. After the 
substitution of Eqs. (|TT j) - (jT?l) and Eqs. (|22 ]) -(|2"3 ]l into Eq. 
(14"T1) . the terms of different orders in N can be identified 
in Eq. fll?}. The leading order terms yield Eqs. ©-([11 
with the substitutions s(t) = s(t),e(t) = e(t),i(t) = i(t): 



ds 

~dt 

de 

~df_ 

di 
dt 



fi(l-s(t)) - (3 (1 + /3 1 cos2irt)s(t)i{t), (25) 
A)(l + Pi cos2irt)s(t)i(t) ~(x + n)e(t), (26) 



X e(t)- (7 + /*)*(*)• 



(27) 



The next-to-leading order terms yield a multivariate lin- 
ear Fokker-Plank equation for the probability distribu- 
tion no,t) HE S3 



d(x k n) 

dxj 



d 2 u 

dxjdxk 



(28) 



where j, k = 1, 2, 3. In Eq. 
Eqs. ©-© 



A(i) 



(-n- mm 

V 



, A(t) is the Jacobian of 

-mm 

-(X + M) /8(*)*(t) ! (2 ( M 
X -(7 + M) 



and B(t) is the symmetric cross correlation matrix 



B(t) 



/ /ii(i) -/12W 
-/«(*) /22W -/23W 

V -/a3(t) /3a(t) 



(30) 



and 



/ u (t) = M(l + «(*)) + /i2(t), 
/aa(t) = /ie(t) + /i 2 (t) + /a3(t), 
/as(«) = (7 + M)i(*) + /2sW 



/ X2 (t) - P(t)s(t)Kt), 
/aa(*) = *e(t). 



(31) 
(32) 
(33) 



(34) 
(35) 



Both matrices A(t) and B(£) are evaluated on the limit 
cycle solutions of Eqs. ©-© and thus they are periodic 
functions of t with the same period T of the limit cycle, 



(36) 



A(t)=A(f + T), B(t) = B(t + T). 



In previous studies of unforced epidemic models [28 
3l| , the power spectrum of stochastic fluctuations Xj (t) 
was computed from the multivariate Langevin equation 
associated with Eq. (|28|) [iolliH 



dt 



Y,A. jk {t)x k {t) + L 3 {t), j,k = 1,2,3, (37) 



where Lj (t) are Gaussian noise terms with zero mean 



and with the correlator 



(L j (t)L k (t'))=B jk (t)6(t-t'). 



(38) 



(39) 



The analytical calculation of the power spectrum through 
the Fourier transform of Eq. (|37|) done in those stud- 
ies depends on the fact that the matrices A and B arc 
constant in the unforced case. For the seasonally forced 
stochastic SEIR model, the matrices A(i) and B(i) are 
time dependent and this method does not work anymore. 
However, one can use the periodicity of A(i) and B(f) 
and Floquet's theory to find a solution of Eq. (|37| and 
compute its power spectrum, as briefly outlined below. 
This method was developed to study the effects of ex- 
ternal noise in nonlinear oscillations close to bifurcations 

43l 44 1 , and it has been applied to intrinsic stochasticity 
for several systems similar to the model considered here 

33]. 

The general solution of the inhomogeneous system of 
first-order linear differential equations (|37|) with matrix 
function A(f ) and vector function L(i) can be written as 
a sum of the general solution of the homogeneous system 



^M = j2A jk (t) Xk (t), 



j,k= 1,2,3, (40) 



and a particular solution of the inhomogeneous system 
45]. Furthermore, the general solution of the homoge- 
neous system (I40[) with A(t) periodic with period T obeys 
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Floquet's theorem 4|| 46 1. This theorem states that if 
X(t) is a fundamental matrix of Eq. (|40"|) . then there 
exists a periodic nonsingular matrix Q(t) with period T 
and a constant matrix R such that 



X(t) = Q(t)e , Q(t) = Q(t + T). 



(41) 



The matrix D = e TR is sometimes referred to as the 
monodromy matrix of the fundamental matrix X(t). Al- 
though D is not unique, its eigenvalues, Ai, A2, A3, called 
the characteristic (Floquet) multipliers associated with 
the periodic matrix A(i), are unique. The eigenval- 
ues of matrix R, pi,p2,P3, are called the characteristic 
(Floquet) exponents associated with the periodic matrix 
A(t). The latter are related to the characteristic multi- 
pliers by 



P] = j, log 



Ail, j — 1, 2, 3, 



(42) 



where the principal value of the logarithm is taken. 

Using further Floquet's theory an analytical expression 
can be obtained for the autocorrelation function of the 
stochastic fluctuations 



Cn(r) 



( Xj {t)x 3 (t + T))dt, j — 1, 2, 3, (43) 



in terms of the matrices Q(t), R and B(t) [33l l43j. Note 
that Eq. (|43|) does not depend on the initial condition 
if the initialization time is set to the infinite past. The 
power spectrum (power spectral density) Pj(cu) of the 
stochastic fluctuations can then be computed from the 
autocorrelation function Cjj (r) via the Fourier transform 
(the Wiener- Khinchin theorem) 



+00 



P» 



Cy(r)e- 



r dr. 



(44) 



In the following section, we will compare the power 
spectra calculated analytically through the method de- 
scribed above with those measured from simulations of 
the stochastic process defined by Eqs. (ITT1) - ([!?)) . The 
spectra of the fluctuations of infectives are of particular 
interest to us because they can be compared with the 
spectra computed from real incidence data 22|, |47 1 . 



IV. RESULTS 

The stochastic process is simulated with the use of a 
modified Gillespie algorithm [48| to account for the ex- 
plicit time dependence in the transition rates 0, l50j |. 
Unless explicitly stated otherwise, all simulations start 
from a random initial condition and time series of 500 
years or 100 years are recorded (50 or 10 years of which 



Pi 


T 


Floquet multipliers 


Floquet exponents 


0.02 


1 


-0.8262 ± i 0.3007 
9.1867 x 10" 60 


-0.1287 ±i 2.7926 
-135.9373 


0.05 


1 


-0.8369 ± i 0.2696 
9.1868 x 10" 60 


-0.1287 ±i 2.8299 
-135.9373 


0.10 


1 


-0.8724 ±i 0.1091 
9.1866 x 10 -60 


-0.1287 ±i 3.0172 
-135.9374 


3 


-0.2715 ±i 0.6119 
7.7134 x 10~ 178 


-0.1338 ±i 0.6628 
-135.9391 


0.12 


2 


0.7618 ±i 0.1310 
8.4390 x 10 -119 


-0.1287 ±i 0.0852 
-135.9374 


3 


-0.6562 ±i 0.1305 
7.7122 x 10" 178 


-0.1340 ± i 0.9878 
-135.9391 


0.2 


2 


0.1184 ±i 0.7630 
8.4357 x 10" 119 


-0.1293 ±i 0.7085 
-135.9376 



TABLE II: Floquet multipliers, Ai,Aa,A3, and Floquet ex- 
ponents, pi,p2,P3, for limits cycles of different periods and 
several values of the forcing amplitude. All other parameters 
are kept fixed as in Table [T] 



are considered as transient). By definition the model does 
not include injection of infectives, thus if an extinction 
occurs before this time a simulation is discarded. The 
largest system size we tested is N = 10 8 and the smallest 
one is N — 5 x 10 5 . For some amplitudes of seasonal forc- 
ing the number of extinctions gets huge and we were not 
able to run such long simulations for system sizes smaller 
than N = 10 6 . However, shorter time series could have 
been obtained (for comparison with the prevaccination 
era these should be about 20 years long jl4|). 

There is one difficulty in the computation of the an- 
alytical autocorrelation functions and analytical power 
spectra. Although an explicit expression for the autocor- 
relation function can be found, the final results have to be 
obtained numerically because the stable limit cycle is not 
known in closed form. In the endemic equilibrium, the 
eigenvalues of the Jacobian of the unforced SEIR deter- 
ministic model ©-(Jll) differ by two orders of magnitude 
for the parameter values typical of childhood infectious 
diseases. This is reflected in the seasonally forced deter- 
ministic SEIR model ©-(31) 17] in the same parameter 
region. For the stable limit cycles given in Table HI two 
of the characteristic multipliers are always complex con- 
jugate and have real part of order 1, while the third one, 
A 3 , is of order 10~ 59 , 10~ 118 and 10~ 177 for limit cycles 



of periods 1, 2 and 3, respectively. In this way, the mon- 
odromy matrix D becomes singular (detD = A1A3A3) in 
double precision numerical calculations and the periodic 
matrix Q(t) cannot be computed. One way to circumvent 
this difficulty is to reduce the computations to the two di- 
mensional central manifold associated with the complex 
eigenvalues, disregarding the dynamics along the strongly 
stable manifold associated with the small real eigenvalue 
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A3 where fluctuations will be strongly damped (a similar 



0.015 



approach was used heuristically in Ref. 22]). Here we 
adopt a direct approach by implementing a Runge-Kutta 
7-8 method for the numerical integration of the differen- 
tial equations on the limit cycle an d by using arbitrary 
precision libraries NTL and GMP 51[. Typically the 



working precision set up in the integrations was 5 dig- 
its higher than the smallest characteristic multiplier (for 
example, 65 digits for limit cycles of period 1) and the 
numerical trajectory was correct up to 50 digits. Thus 
we could perform all the computations in the original 
variables of the full three dimensional system, which al- 
lows an immediate comparison with the power spectra of 
the simulations. The summary of the Floquet multipliers 
and Floquet exponents computed in this way is given in 
Table HI 

Note that until now we have not raised the question of 
the stability of the deterministic limit cycles because we 
do know from previous studies that the cycles are stable 
[tR fl9l l2p| . However, we automatically check this by 
computing the absolute values of Xj which must be less 
than unity for a limit cycle to be asymptotically stable 
(or , alternatively, the real parts of pj must be negative) 
[45|,[46j]. The homogeneous equation (|40|) is, in fact, a 
variational equation for small perturbations Xj (t) around 



the periodic solution of Eqs. ©-Q [45j, |46j. The pres- 
ence of two complex conjugate Floquet multipliers im- 
plies that deterministic perturbations decay to the cycle 
in a damped oscillatory way. This situation is similar 
to that of unforced systems with a stable focus in which 
resonant amplification of stochastic fluctuations was ob- 
served 
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291 ] . It will become clear in what follows that 
in the stochastic system with seasonal forcing fluctua- 
tions around the limit cycles are significantly amplified 
too, and that the analytical and simulated power spectra 
cannot be explained by the deterministic theory alone. 

To validate the theory developed in this study, we first 
compare analytical and simulated autocorrelation func- 
tions for the stochastic fluctuations of infectives. Typical 
values of the forcing amplitude are chosen in the annual 
and biennial regime (see left and right upper plots of Fig. 
[2|). Far from bifurcation points and for large system sizes, 
we find an excellent agreement in both cases (a small di- 
vergence can occur because of the sparse discretization 
of the orbit which we are lead to do to avoid very heavy 
computations). The lower plots of Fig. [2] show the corre- 
sponding power spectra. More exactly, the x-axes stands 
for temporal frequency and the y-axes is the power of the 
discrete Fourier transform of the autocorrelation func- 
tion time series. In this and in the following plots, the 
power spectra are normalized so that the total power is 
the summed squared amplitude of the time series [52| . 
One observes that the power spectra exhibit a number 
of peaks occurring regularly. As it has been noticed be- 
fore by several authors [33, ,43[, the peaks are located at 




frequency (year ) 



frequency (year ) 



FIG. 2: (Color online) The upper panels show the autocor- 
relation functions of the stochastic fluctuations of infectives. 
The theoretical curves are plotted in gray (blue), and the 
curved computed from the simulations are plotted in black 
(for the annual cycle the analytical and simulated autocor- 
relation functions are strictly coincident). The lower panels 
show the corresponding power spectra. The vertical helper 
lines mark the frequencies predicted by Eq. (|45[) . The dashes 
(dotted) lines are calculated by taking the plus (minus) in the 
equation. The system size used for simulation is 10 8 . 



frequencies 



n |Im pi ;2 | 
T 2tt 



(45) 



where n = 0, ±1,±2, |Im pi^\ denotes the absolute 
value of the imaginary part of complex conjugate Floquet 
exponents (see Table HB. and T is period of a limit cycle. 
Note that for the annual limit cycle the main stochastic 
peak is situated at |Im jc>i,2 1/ (27r) while for the biennial 
cycle this peak is much smaller than the peak at 1/2 — 
|Im pi,2|/(27r) which is dominant. The plots are done in 
lin-log scale for a better observation of the structure of 
the power spectra. 

The analytical predictions for stochastic power spectra 
work well in a large range of /3i for annual and biennial 
limit cycles. However, as the system size is decreased sys- 
tematic deviations start to appear. We first consider the 
case of the annual limit cycle induced by a low forcing am- 
plitude. In Fig. [3] we compare the power spectra for the 
number of infective individuals from simulations for sev- 
eral system sizes N. In all panels, sharp peaks due to the 
annual limit cycle can be observed at integer frequencies, 
as well as several broader stochastic peaks whose frequen- 
cies are given by Eq. (|4"B"]) . For all system sizes the first 
(and largest) stochastic peak is situated at |Im pi,2 |/(27r), 
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frequency (year ) frequency (year ) 



FIG. 3: Power spectra of the number of infectives calculated 
from simulations for several system sizes N. The simulations 
of 500 years were run for N = 10 s , 10 7 , 5 x 10 6 , 10 6 and of 100 
years for N = 5 x 10 5 . Averages over 10 3 realizations were 
made to obtain each curve. The forcing amplitude is low and 
corresponds to the annual limit cycle, /3i = 0.02. The vertical 
helper lines mark the frequencies predicted by Eq. (|45|) . 

the second is situated at 1 — |Im pi,2 1 / (27r) , etc. However, 
for N — 10 6 and N = 5 x 10 5 the dominant frequencies of 
the stochastic peaks are very slightly shifted to the left, 
so the characteristic period of the stochastic fluctuations 
gets higher. This effect has already been observed before 
in the spectra of fluctuations around a stable node, in 
unforced systems [28j]. We also find that the spectra arc 
fully described by the theory only for very large popula- 
tions. Beginning with N — 10 7 a number of much smaller 
stochastic peaks starts to appear close to the determin- 
istic peaks. Their amplitudes increase as the population 
size decreases but they stay orders of magnitudes lower 
than the dominant stochastic and deterministic peaks. 
These secondary peaks cannot be explained within the 
theory developed in this study and require considering 
corrections to the linear Fokker-Planck equation. An- 
other effect apparent in small systems is the change in 
the relative amplitude of the main stochastic and deter- 
ministic peaks. For small populations the deterministic 
limit cycle does not dominate the dynamics of the sys- 
tem any more. The enhancement and broadening of the 
stochastic peaks show a much noisier and irregular dy- 
namics. For comparison, see the panel for N = 5 x 10 D 
in Fig. |3] where the main stochastic peak is significantly 



FIG. 4: Power spectra of the number of infectives calculated 
from simulations for several system sizes N. The simulations 
of 500 years were run for N = 10 8 , 10 7 , 5 x 10 6 , 10 6 and of 100 
years for N — 5 x 10 5 . Averages over 10 4 (10 3 ) realizations 
were made to obtain the power spectra for N = 10 8 (for the 
remaining system sizes). The forcing amplitude is 2.5 higher 
than in Fig. [3] but the deterministic system is still in the 
annual limit cycle regime, /3i = 0.05. The vertical helper 
lines mark the frequencies predicted by Eq. (|45[) . 

higher and broader than the main deterministic peak. 

The same picture of changes in the spectra for different 
population sizes is maintained if the amplitude of forcing 
is increased but the annual limit cycle is still stable. As 
an example, see Fig. 2] done for fix = 0.05 (this forcing 
amplitude is 2.5 larger than the one used in Fig. [3]). 

If the forcing amplitude is increased even further the 
period doubling of the limit cycle occurs at (3[ « 0.11479 
and everywhere in the interval 0.11479 < /?i < 0.2 a sta- 
ble limit cycle of period 2 is present in the deterministic 
model. In Fig. [5] we compare the spectra from simu- 
lations and the analytical spectra for a range of system 
sizes. Again, the spectra demonstrate narrow determin- 
istic peaks at multiples of 1/2 due to the limit cycle of 
period 2 and regular stochastic peaks. For the largest 
system size, the stochastic peaks are predicted correctly 
by the theory. The major stochastic peak is now lo- 
cated at 1/2 — |Im pi,2 1/ (27r). Note that the peak which 
used to be dominant in the annual regime was located 
at |Im pi,2|/(27r). But now this frequency corresponds 
to the smallest among all stochastic peaks in the range 
of frequencies shown in the plot. This observation de- 
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N=10' 



frequency (year ) 



frequency (year ) 



FIG. 5: Power spectra of the number of infectives calculated 
from simulations for several system sizes N. The simulations 
of 500 years were run for N = 10 8 , 10 7 , 5 x 10 6 and of 100 
years for N = 10 6 . Averages over 5 x 10 3 (10 3 ) realizations 
were made to obtain the power spectra for N = 10 s (for the 
remaining system sizes). The forcing amplitude corresponds 
to the biennial limit cycle, /3i = 0.2. The vertical helper lines 
mark the frequencies predicted by Eq. (|45l) . 



serves a longer comment and we will go back to it later. 
For smaller values of N we observe that the determin- 
istic peak at 1/2 gets smaller and merges with the two 
surrounding stochastic peaks. At some point it becomes 
impossible to distinguish between the stochastic and de- 
terministic peaks which are represented by a single broad 
peak around 1/2. The stochastic peaks get significantly 
enhanced but also the range of frequencies present in the 
infective time series. 

The next parameter range of interest is close to the 
point where the period doubling bifurcation occurs for 
the deterministic system. In this regime we do not expect 
to have an agreement between the theory and simulations 
because none of the cycles is stable enough, however there 
are interesting conclusions which can be drawn from the 
power spectra. Recall that the period doubling in the 
deterministic model occurs at f3[ 0.11479 but in the 
stochastic model the transition from the annual limit cy- 
cle to the biennial limit cycle is shifted to a higher value 
of In Fig. [6] we compare the spectra for the number of 
infectives from simulations performed for two close values 
of P\ in the vicinity of the deterministic period doubling 
point, one before the bifurcation (/3i = 0.1, correspond- 
ing to the stable annual limit cycle, gray (blue) curves) 
and one after (/3i = 0.12, corresponding to the stable 
biennial limit cycle, black curves). As one can see the 
transition is blurred and shows up later in the stochas- 
tic system. The same behavior of the fluctuations power 
spectrum around a bifurcation was described in Ref. [43j j . 




0.5 1 1.5 2 2.5 
frequency (year 1 ) 



FIG. 6: (Color online) Power spectra of the number of in- 
fectives calculated from simulations for several system sizes 
N. The gray (blue) and black curves were obtained for 
/3i =0.1 and /3i = 0.12, respectively. In the determinis- 
tic model the period doubling occurs at /3j ~ 0.11479. The 
frequencies of the vertical helper lines correspond to the an- 
nual limit cycle. The simulations of 500 years were used for 
N = 10 8 , 10 7 , 5 x 10 6 , 10 6 and of 100 years for N = 5 x 10 5 . 
Averages over 10 3 realizations were made to obtain all curves. 
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n=io/ 



frequency (year ) 




frequency (year ) 



FIG. 7: (Color online) Power spectra of the number of infec- 
tives calculated from simulations for two system sizes N and 
two sets of initial conditions [for the gray (blue) curve sim- 
ulations started from random initial conditions and for the 
black curves the initial conditions were chosen close to the 
deterministic triennial cycle]. In the left (right) panel the 
vertical helper lines mark the predicted peak frequencies for 
the triennial (annual) limit cycle. The amplitude of seasonal 
forcing corresponds to coexisting stable annual and triennial 
limit cycles in the deterministic model, /3i = 0.1. 
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There is almost no difference between the simulated spec- 
tra. The first deterministic peak signaling the biennial 
limit cycle has to be present at 1/2 for the black curves, 
and this is where the stochastic peaks for the annual limit 
cycle are located. The two main stochastic peaks of the 
annual cycle move on to 1/2 and at some point give rise 
a deterministic peak of the biennial cycle. This is seen 
from the plot for the largest system size, however the de- 
terministic peak becomes obvious if the pi is increased 
a bit further. A similar broadening and shifting of the 
bifurcation point in simulations has also been observed 
in the unforced model with a transition from a stable 
focus to a stable limit cycle through an Andronov-Hopf 
bifurcation [ji.!^. 

The last question one naturally poses is how does the 
coexistence of attractors influence the power spectrum, 
through possible switches between basins of attraction 
induced by stochasticity. The relevance of this mecha- 
nism has been argued for in the literature 



53, |54| 



particular to try to provide an explanation for the di- 

" SH- In the 
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versity of patterns of childhood diseases 
deterministic model there is a range of j3i where stable 
annual and biennial limit cycles coexist with a stable tri- 
ennial cycle. In particular, this happens for the parame- 
ter values of Fig. [6] We, however, have not identified any 
deterministic peaks at multiples of 1 /3 associated with a 
triennial limit cycle or even an indication of such a be- 
havior in the stochastic system. The same was confirmed 
for other values of /3\ where the orbits coexist. To clar- 
ify this situation we performed an experiment in which 
one set of simulations started on a triennial limit cycle 
(these conditions are favorable for the triennial cycle as 
it will become clear later) and the other set started from 
a random initial condition. The resulting power spectra 
are shown in Fig. [7] for /3% — 0.1. The black [respectively, 
gray (blue)] curves are the power spectra of the infec- 
tive time series beginning from favorable (respectively, 
random) conditions. It is only for system size N = 10 8 
that we were able to observe a triennial limit cycle with 
the fluctuations around it described by the van Kampen 
expansion (see the left panel in Fig. the helper lines 
mark the frequencies for the triennial limit cycle). For 
all other system sizes the power spectra are identical to 
those already shown in Fig. [5] 

In this case, a better insight into what happens in the 
simulations is given by inspecting the time series. In Fig. 
E] we show the density of infectives in a typical single 
realization of the stochastic model starting from the tri- 
ennial limit cycle (for the same value of /?i as in Fig. 
[7]). For system size N — 10 s the density exhibits reg- 
ular triennial oscillations. Their characteristic features 
are a high amplitude (at least twice as large as that of a 
typical annual or biennial cycle) and very low number of 
infectives between recurrent epidemics. For N = 10 s the 
system tends to stay on the triennial limit cycle because 
the fluctuations are not large enough to drive it to the 




FIG. 8: The infective density recorded from a typical realiza- 
tion of the stochastic model starting on the deterministic tri- 
ennial limit cycle. The seasonal forcing amplitude is /3i = 0.1. 
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FIG. 9: The infective density recorded from a typical real- 
ization of the stochastic model starting on the deterministic 
annual limit cycle. The seasonal forcing amplitude is /?i = 0.1. 
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n 

pi 


N = 10 s 


TV = 10 7 


5 x 10 < JV < 10 


0.05 


well defined dominant annual 
deterministic peak 


well defined dominant annual 
deterministic peak and two 
subdominant close to biennial 
broad stochastic peaks 


well defined annual determinis- 
tic peak and close to biennial 
broad stochastic peak 


0.2 


well defined dominant annual 
and biennial deterministic 
peaks 


well defined dominant annual 
deterministic peak and subdom- 
inant broad biennial stochastic 
peak 


well defined dominant annual 
deterministic peak and subdom- 
inant close to biennial much 
broader stochastic peak 



TABLE III: Qualitative features of the power spectra for different values of Pi and TV. 



annual limit cycle. But for smaller TV the system either 
goes extinct or drops very quickly to the annual limit cy- 
cle (or a biennial one if /3i > 0.11479) and stays on it. 
For the population of one million of individuals not even 
one oscillation over the triennial cycle is followed by the 
stochastic system because the relative fluctuations are 
large at the first drop of infectives and drive the system 
to the basin of attraction of the annual cycle. In Fig. [9] 
we show some of the time series for the same parameter 
values as in the previous figure but for initial conditions 
starting on the annual limit cycle. The density exhibits 
alternating regions of annual and biennial oscillations for 
large population sizes with an ever increasing role of the 
stochastic fluctuations as TV decreases. For small sizes 
the time series look quite irregular, occasionally exhibit- 
ing interepidemic intervals longer than 2 years. This case 
corresponds to the power spectra where the stochastic 
peaks are broad and high. It is remarkable that we have 
never observed the backward switching from an annual 
(or biennial) limit cycle to a triennial one. 

The qualitative features of the power spectra for two 
different values of the forcing amplitude, ft\ , and different 
population sizes, TV, are summarized in Table ITTTl As dis- 
cussed above, well defined high amplitude peaks show up 
in the time series as regular cycles, and broad stochastic 
peaks as lower amplitude noisy oscillations. 



V. DISCUSSION AND CONCLUSIONS 

In conclusion, in this study we have considered the 
deterministic and the stochastic SEIR models with sinu- 
soidally varying contact rate. Depending on the forcing 
amplitude, the deterministic model exhibits a range of 
stable attractors, the most visible of which are limit cy- 
cles of periods 1, 2 and 3. Using the van Kampen's expan- 
sion of the master equation of the corresponding stochas- 
tic model we have calculated autocorrelation functions 
and power spectra of the stochastic fluctuations around 
these attractors. We have compared the analytical re- 
sults with those obtained from direct simulations of the 
stochastic model. We have found that in a large range of 
the forcing amplitude there is an excellent agreement be- 
tween the theory and simulations. The prerequisites for 



this are large system sizes (typically higher than 10 6 ) and 
the stability of attractors, namely the more stable the 
cycle and the higher the system size are the better the 
agreement between the simulated and analytical power 
spectra is. This is exactly what we would expect for the 
quality of the approximation given by the truncation at 
first order of the full van Kampen expansion. 

The power spectra of the infective time series demon- 
strate peaks of two types. The narrow peaks are due to 
a limit cycle of a given period and the broader peaks are 
due to the resonant amplification of stochastic fluctua- 
tions around the limit cycle. It has been argued that 
the presence and position of deterministic and stochas- 
tic peaks in the power spectra obtained from real data 
records of childhood diseases can be predicted by the 
deterministic theory alone, and that, moreover, the fre- 
quency of the stochastic peak is defined by the frequ ency 



of the transients near a stable limit cycle [22 



47] 



We have identified the main frequencies of the stochas- 
tic peaks in the annual and biennial regimes and shown 
that these do not necessarily equal the frequency of 
the damped oscillations of deterministic perturbations 
around a cycle. Thus, neither the full structure of the 
power spectrum of the infective time series nor the most 
prominent frequency of the stochastic peak can be fully 
predicted by the deterministic model. 

Another conclusion concerns the role of the coexis- 
tence of attractors in the seasonally forced stochastic 
and deterministic SEIR models. The coexistent stable 
limit cycles present in the deterministic epidemic mod- 
els have been conjectured to be the reason why irregu- 
lar dynamics is observed in the corresponding stochastic 
system. Namely, it has been a systematic assumption 
of severalpapers on childhood infectious diseases mod- 
eling [l9L Ell . |22| based exclusively on the deterministic 
analysis that for small populations the stochastic system 
should constantly switch between cycles of different pe- 
riods due to demographic or environmental noise, thus 
exhibiting irregular dynamics. The results of the present 
study contradict this view. For values of fi\ close to the 
period doubling bifurcation value, demographic as well as 
environmental noise can promote switching between pe- 
riod 1 and period 2 cycles. However, for the stable limit 
cycles of period 1 or 2 coexisting with a stable cycle of pe- 
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riod 3 such switching is absent in the stochastic system. 
Notwithstanding the basin of attraction of the triennial 
limit cycle being roughly 25% of the total initial condi- 
tions [19( for the deterministic system, it is extremely 
unstable in the simulations. Even for large system sizes 
the system is brought to a cycle of lower period after a 
very short time, and backwards switching is not observed. 
This can be understood from the shape of the periodic 
solutions of the triennial cycle of the deterministic model, 
together with the fact that most initial conditions with 
very low number of infectives are attracted to the lower 
period cycles. 

We have become aware of a recent paper where the fluc- 
tuation power spectrum of a seasonally forced stochastic 
SIR model is computed through the same method [55j |. 
In this paper, a transition representing injection of in- 
fectives is included in addition to the processes of infec- 
tion, recovery, and birth-death. It is shown that, even 
for very low infective immigration rates, the inclusion of 
this term has a drastic impact on the bifurcation diagram 
of the deterministic model, leaving stable annual and bi- 
ennial limit cycles as the only possible attractors [HBj ]. 
This shows that the competing higher period attractors 
of the forced system are very fragile indeed. They are 
not robust, in the deterministic setting, with regard to 
small changes in the model. On the other hand, in the 
stochastic setting, their effective basin of attraction is 
negligibly small, except for unrealistically large popula- 



tion sizes. The results of this study and those of Ref. 55 1 
concur to support the view that the main ingredient to 
understand the observed incidence patterns of childhood 
infectious diseases is stochastic amplification, rather than 
noise induced switching between competing attractors of 
the deterministic system. 
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APPENDIX 

As described in the Section III of the main text, the 
master equation for the seasonally forced stochastic SEIR 
model is constructed by considering all possible transi- 
tions increasing or decreasing the probability of finding 
a system in state with mi susceptible, mi exposed and 
7713 infective individuals. This equation reads as follows: 



dV(m 1} m2,m 3 ,t) 
dt 



7X+™mT-l,m 3 ^( m l + 1, TO 2 " 1, m 3 , t) + 7X I ,C^+l' P ( m l> m 2> TO 3 + 1> f ) 

V?um2+7?m 3 -l' P ( m l> m 2 + 1, m 3 - 1, *) + 'C l 1 -T,mT,m 3 ' P ( m l ~ m 2> TO 3, *) 
< C 1 +T,m™m 3 :P ( m l + !' m 2> m 3> f ) + Vni)m2+l,m 3 ^( m l) m 2 + *i m 3, t) 



/mi , m 2 ,' m 3 



T 



s-rm-i ,m 2 , m 3 — 1 
'mi, m 2 ,?7l3 



r 



.. : / , ':</M-"'2-'": i + M) 

+ 

V(mi,m 2: m 3 ,t) 



rym i — 1 , fH2 + 1 , mz 



+ r r . 



mi ,?7l2 — 1,7713+1 



mi ,7772 ,7773 



'■ymi + l, 7772,7773 i <T~?77i —1,7772 ,7773 
' 7771 ,7772 ,7773 ' 7771 ,7772 ,7773 



-y-777 1,777 2 — 1,7773 
' 7771 ,7772 ,7773 



y777i ,7772 ,7773 — 1 



7771 ,7772 ,777 3 



r 



(46) 



The above equation and its subsequent analysis can in the main text by Eqs. (|19|) -(|21 j) . The substitution of 
be simplified with the use of the step operators defined Eqs. ([T9 |) -pT |l into Eq. (|46)) gives 



dV{mx, m 2 , m 3 , t) 
dt 



( — 1 _ -|\ <T-mi,m 2 -l,m3 + l 
1^3 / 'mi,m2,m3 



— 1,7772,7773 
7772,7773 



+ (e 3 -l)-C- 



7772 ,7713 — 1 
7772,7773 



/ -1 _ l \ q-mi- 1,7772 + 1,7773 , / - 1 _ 1 \ ^7771+1,7772 , 7773 , / _ -1 \ ^777 1 , 777 2 - 1 , 777 3 

V 1 2 / ' 7771,7772,7773 ' V I / 7771,7772,7773 ~ \ 2 ^ ) ' 7771,7772,7773 

V(mi,m2,m 3 ,t). 



I ^) ' 1711,012,013 



(47) 



Press, Princeton, 2007). 
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